/*
 * test.cpp
 *
 *  Created on: Nov 4, 2010
 *      Author: jagiella
 */

#include <stdio.h>

#include "SparseMatrix.h"
#include "Solver.h"


int main( int argc, char **argv)
{

	int n=1001;

	Solver<double> *S = new Solver<double>( n, Solver<double>::BiCGSTAB);
	S->maxit = 1000;
	S->maxerr= 1e-20;

	SparseMatrix<double> *A = new SparseMatrix<double>(n,n);
	SparseMatrix<double> *B = new SparseMatrix<double>(n,n);
	double b[n];
	double x[n];

	// INIT SYSTEM
	fprintf( stderr, "INIT SYSTEM\n");
	for( int i=0; i<n; i++) x[i]=0;
	(*B)(0,0,1);(*A)(0,0,1); b[0] = 0;
	for( int i=1; i<n-1; i++){
		(*B)(i,i-1,-1);(*A)(i,i-1,-1);
		(*B)(i,i  , 2);(*A)(i,i  , 2); b[i] = 1;
		(*B)(i,i+1,-1);(*A)(i,i+1,-1);
	}
	(*B)(n-1,n-1,1);(*A)(n-1,n-1,1); b[n-1] = 0;
	/*A->printMatrix( "M", "%10.3lf");
	fprintf( stderr, "b = [\n");
	for( int i=0; i<n; i++){
		fprintf( stderr, "%10.3lf", b[i]);
	}
	fprintf( stderr, "\n]\n");
*/

	S->IncompleteLUfactorization(B);
	double temp[n];
	S->SolveLowerTriangular( B, b, temp);
	S->SolveUpperTriangular( B, temp, x);


	// SOLVE SYSTEM
	fprintf( stderr, "SOLVE SYSTEM\n");
	S->Solve(A,b,x);
	/*fprintf( stderr, "x = [\n");
	for( int i=0; i<n; i++){
		fprintf( stderr, "%10.3lf", x[i]);
		x[i]=0;
	}
	fprintf( stderr, "\n]\n");*/

	fprintf( stderr, "LU DECOMPOSITION\n");

	S->IncompleteLUfactorization(A);

	//A->printMatrix( "M (LU)", "%10.3lf");

	double y[n];
	double xILU[n];

	fprintf( stderr, "SOLVE SYSTEM: Ly=b\n");
	//S->SolveLowerTriangular( A, b, y);

	fprintf( stderr, "SOLVE SYSTEM: Ux=y\n");
	//S->SolveUpperTriangular( A, y, xILU);

	S->SolveLU( A, b, y, xILU);

	for( int i=0; i<n; i++)
		fprintf( stdout, "%i %lf %.7lf\n", i, x[i], xILU[i]);

	/*fprintf( stderr, "y = [\n");
	for( int i=0; i<n; i++){
		fprintf( stderr, "%10.3lf", y[i]);
	}
	fprintf( stderr, "\n]\n");

	fprintf( stderr, "x = [\n");
	for( int i=0; i<n; i++){
		fprintf( stderr, "%10.3lf", x[i]);
		x[i]=0;
	}
	fprintf( stderr, "\n]\n");
	 */

	return 1;

	/*double x[5];
	double b[5];

	A->set(0,0, 1);	b[0] = 1;
	for(int i=1; i<4; i++){
		A->set(i,i-1, -1);
		A->set(i,i  ,  2);
		A->set(i,i+1, -1);
		b[i] = 1;
	}
	A->set(4,4, 1);	b[4] = 1;

	Solver<double> *S = new Solver<double>( 5, Solver<double>::CG);
	S->maxit = 10000;
	S->maxerr= 1e-20;

	int it = S->Solve( A, b, x);

	fprintf( stderr, "\nSolution after %i iterations:\n", it);
	for(int i=0; i<5; i++)
		fprintf( stderr, "%lf\n", x[i]);

	return 0;*/
}
